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Phone: (216) 433-8334 
Fax: (216) 433-8643 
Email: dpaxson@grc.nasa.gov 


Abstract 

A one-dimensional, CFD based combustor simulation 
has been developed that exhibits self-excited, thermo- 
acoustic oscillations in premixed combustor geometries 
that typically have large, abrupt changes in cross 
sectional area. The combustor geometry is 
approximated by dividing it into a finite number of one- 
dimensional sectors. Within each sector, the equations 
of motion are integrated numerically, along with a 
species transport and a reaction equation. Across the 
sectors, mass and energy are conserved, and momentum 
loss is prescribed using appropriately compatible 
boundary conditions that account for the area change. 
The resulting simulation and associated boundary 
conditions essentially represent a one-dimensional, 
multi-block technique. Details of the simulation code 
are presented herein. Results are then shown comparing 
experimentally observed and simulated operation of a 
particular combustor rig that exhibited different 
instabilities at different operating points. It will be 
shown that the simulation closely matched the rig data 
in oscillation amplitudes, frequencies, and operating 
points at which the instabilities occurred. Finally, 
advantages and limitations of the simulation technique 
are discussed. 

Introduction 

The problem of thermo-acoustic combustion 
instabilities in aircraft and land-based gas turbine 
engines is currently the focus of numerous research 
efforts 1 . Lean-premix, low emission combustor designs 
appear to be particularly susceptible to this 
phenomenon due to, among other things, their lack of 
cooling air and their acoustically stiff geometry. One 


+ Senior Member, AIAA. 

Copyright © 2000 by the American Institute of Aeronautics and 
Astronautics, Inc. No copyright is asserted in the United States under 
Title 17, U.S. Code. The U. S. Government has a royalty-free license to 
exercise all rights under the copyright claimed herein for Governmental 
Purposes. All other rights are reserved by the copyright owner. 


possible solution to the instability problem is active or 
feedback control. Here, a signal from the system (e.g. 
pressure) is used to detect incipient instabilities, and 
some form of actuation (such as fuel flow perturbation) 
is used to suppress and restabilize the system. This is in 
contrast to present passive methods where the structure 
of the combustor is physically changed to avoid the 
instability entirely. Successful active control design 
however, is greatly enhanced by accurate modeling and 
simulation of the combustor of interest. The essential 
physical phenomena should be correctly captured. For 
example, the simulation should self-excite, as is often 
the case with the actual instability, as opposed to 
requiring some form of external forcing. On the other 
hand, dynamic characterization of the system required 
for control design necessitates parametric analyses and 
multiple simulation runs. This, in turn, places the 
practical requirement of high speed on the simulation. 

For many combustor configurations, the underlying 
physics of the system are simply too complex to lend 
themselves to any sort of simplified approach that is 
practical for active control design. Instabilities that are 
both longitudinal and circumferential for instance, or 
those that are intricately tied to some vortex shedding 
phenomenon are examples. 

If however, the instabilities are predominantly oriented 
along one spatial dimension, and if consideration is 
limited to premix type combustion, then it may be 
possible to satisfy the simulation objectives of 
numerical speed (i.e. model simplicity), and physical 
accuracy. Such restrictions, though severe, apply to a 
large number of combustors. This paper introduces a 
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Sector 1 

Sector 2 


Injector Region 


Figure 1. — Lean premix combustor schematic. 

method to achieve the simulation goals just described 
with the added, essential, capability of accounting for the 
large and abrupt changes in cross sectional area typically 
associated with premix combustors. An example of such 
geometry is shown schematically in Fig. 1 . 



Sector 3 


Combustor Region 


frequencies, amplitudes, and operating points. The 
paper concludes with a discussion of the possible 
advantages and limitations of the simulation technique. 

Governing Equations 

Given a combustor geometry which has been suitably 
divided into one-dimensional sectors, within each sector 
the numerical integration is performed on the following 
governing differential equations for a caloric ally perfect 
gas, written below in non-dimensional form 


3w 3F(w) 
dt 3x 


S(w,x) 


where 


( 1 ) 


The combustor geometry is approximated by dividing it 
into a finite number of one-dimensional (constant area) 
sectors. Within each sector, the equations of motion are 
integrated numerically, along with a species transport 
and a very simple one-step reaction equation. This 
aspect of the simulation has been used successfully to 
investigate both deflagrative and detonative combustion 
processes for a variety of applications 2 4 . Across the 
sectors, mass and energy are conserved, and momentum 
loss is prescribed using appropriately compatible 
boundary conditions that account for the area change. 
This approach avoids the difficulties often found when 
quasi-one-dimensional computational schemes are 
applied to geometries with abrupt area changes*. 
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Details of the simulation methodology are presented 
herein. These include a presentation of the assumed 
governing equation set, a description of the numerical 
integration scheme, and an algorithmic description of 
the boundary conditions that tie two sectors together. A 
simple, non-reactive acoustic flowfield simulation is 
then shown which serves (somewhat heuristically) to 
validate the sectored approach. A Pulse Jet simulation is 
then presented which clearly demonstrates the 
capability of capturing thermo- acoustic oscillations in a 
device that exhibits them by design. Comparisons are 
then made between experimentally observed and 
simulated instabilities in a specific combustor rig. The 
rig was a critical test of the approach in that it had 
multiple, abrupt area changes and it exhibited two 
different unstable limit cycles at two different operating 
points. It will be shown that the simulation successfully 
exhibited the self-excited oscillations at the correct 


The distance, x has been normalized by the combustor 
length, L. The time, t has been normalized by the 
characteristic wave transit time, Ua \ where a* is the 
speed of sound at a chosen reference state. The 
pressure, p and density, p have been normalized by their 
respective reference values and the axial velocity, u has 
been normalized by a . The mass fraction of reactant is 
Z . Note that z can be related to, but is not the same as, 
actual fuel fraction. The ratio of specific heats is 
denoted by y 

For this non-dimensionalized form of the equations, the 
equation of state is written 

P = PT (4) 

The speed of sound is simply Vt. 


* These arc discussed briefly in the appendix of this paper. As of this 
printing, the author has not found a q-l-d approach that successfully 
simulates the particular combustor presented herein. 
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The non-dimensional source vector is written as: 


S(w,x) = 


It may contain contributions from reaction, turbulent eddy 
diffusion, wall viscous forces, heat transfer, and 
supplementary flows (including reactant flow such as a 
fuel injector). The Reynolds number, Re is defined as 
p a Up. The turbulent viscosity ratio, § is defined as the 
ratio of turbulent to molecular viscosity, p t /p. The heat of 
reaction is q 0 , which has been normalized by the square of 
the reference speed of sound. R is the reaction rate, 
defined below. Pr t and Sc t are the turbulent Prandtl and 
Schmidt numbers respectively. It is noted that the ratio p { 
//i, as well as Pr t and Sc f , are not necessarily physically 
based. They are chosen to represent 'effective' diffusion 
coefficients which allow the diffusion terms of the 
governing equations to mimic the transport behavior of, 
say, large recirculation zones and/or shear layers. For the 
results presented in this paper, the value of the wall 
friction coefficient, C 2 is zero. 
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Combustion Model and Reaction Rate 

The reaction rate of the modeled combustion process is 

a very simple, single species equation of the form: 


R = K 0 pz(C, - C 2 z> 


l-fen/T.) ;Tj > T ig 


•T < T 

* A i A ig 


( 6 ) 


where Ko is the reaction rate constant (normalized by 
the reference wave transit time), and t^ 2 are constants 
defining the type of reaction, T igI1 is the ignition 
temperature and Tj is the temperature in the i* 
numerical cell. This reaction model has been used with 
success for, among others, the numerical simulation of 
propulsion devices with non-steady burning 2 . For the 
results presented herein, the value of C, 2 was set to 0.0. 


Flameholding 

A model for flameholding is not shown per se in the 
governing equations. It is modeled as a distribution in 
the value of the turbulent diffusion coefficient, £ t . In 


particular, the value in the simulation results to be 
presented changes from zero to some maximum 
prescribed value over 5% of the combustor length. It is 
then gradually reduced to zero over the distance 
between the location of the maximum value and the end 
of the sector. The location of this transition must also be 
prescribed. Typically, it begins in the vicinity of a 
known dominant recirculation zone. Obviously, this is 
an extremely crude model for flameholding and 
currently requires a degree of 'tuning' to avoid such 
phenomena as blow-out. Furthermore, the particular 
form and location of the diffusion coefficient 
distribution has been found to strongly influence the 
formation of self-excited oscillations. On the other 
hand, the notion of using a diffusion equation to model 
recirculation and mixing is a natural one that, with 
additional analysis, may be possible to generalize and 
quantify. 

Heat Transfer 

The term Q ht in Eqn. 5 is represented by a simple 
algebraic expression 

Qm = “(Tjnf -Tj) (7) 

where a is a user specified constant (and may depend 
on the velocity and density if convective), T inf represent 
an assumed temperature of the surface or space to or 
from which heat may be transferred, and T* is the 
(computed) gas temperature in the i A cell. In the 
combustor simulation results to be presented, this term 
is used to model the effects of cooling spray in the 
downstream portion of an experimental rig. The value 
of a, and the combustor extent over which the term acts 
was determined by matching simulation results with 
experimental measurements. 

Numerical Method 

The simulation numerically integrates the above 
equations of motion using a very simple, second- order 
MacCormack scheme 5 . A Baldwin-MacCormack type 
artificial viscosity scheme has been added in order to 
damp non-physical oscillations in vicinities of strong 
spatial gradients such as those brought about by the 
combustion process. Typically, this type of artificial 
viscosity uses a spatial pressure profile as a weighting 
coefficient. In this simulation a density profile is used 
instead for the mass, momentum, and energy equations, 
and a reaction fraction profile is used in the species 
transport equation. The numerical scheme and 
associated artificial viscosity were chosen because of 
the computational speed which they afford. 
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Boundary Conditions 

Boundary conditions are the means by which area 
changes are accounted for, and by which numerical 
image cell states are prescribed. In hyperbolic systems 
they allow computation of interior cell states. Boundary 
conditions at each sector may be imposed as either 
partially opened, fully open, choked inflow (i.e. 
constant mass flux), or solid walls. In any case the code 
anticipates the flow direction and applies appropriate 
(i.e. well-posed) states to image cells just outside the 
computing domain. These are described briefly below. 
Details may be found in Refs. 6 and 7. 

Fully Open 

If the flow is outward from the computing domain, only 
the static pressure is imposed. The remaining 
information (density, velocity, and reaction fraction) is 
obtained via characteristic jumps across the adjacent 
numerical cell. If the flow is inward, total pressure and 
temperature, and reaction fraction are imposed. The 
remaining unknown quantity (velocity) is obtained 
through iteration and characteristic jumps across the 
adjacent numerical cell. In this case the image cell state 
ultimately meets the simultaneous requirements of 
having the imposed total pressure and temperature, and 
having no outward-running characteristics. 


Partially Open 

This type boundary is depicted in Fig. 2. Depending 
upon the external conditions, the flow will be either in or 
out of the computing domain. Looking first at the inflow 
scenario (a), the external flow is modeled as isentropic 
and steady relative to some prescribed stagnation state (as 
in the fully-open case above). Adjacent to this region is 
an imaginary "mixing region" where the fluid changes 
from a highly non-uniform distribution across the 
partially open tube on the left to some mixed out average 
on the right. The mixing is modeled as instantaneous, and 
the extent of the mixing region is modeled as 
infinitesimal. Adjacent to this region is the first "cell" of 
the computing domain for which, at the present time, the 
conditions are completely known. The two zones are 
related by the following mass momentum and energy 
equations: 


p 0 uo = P e u e — 

At 


Po =P e -Y 


PeU e 


l 

A, 


fjr ' 

U-ij 


UoPo + ~Po uo 3 -Pe u e H — 
1 A, 


( 8 ) 

(9) 

( 10 ) 


Here, H is the total enthalpy. The subscripts 0 and e in 
these equations designate the end and the beginning 



e-plane 

(a) 

isentropic zone 



(b) 

Figure 2. — Partially open boundary condition 
schematic; (a) Inflow, (b) Outflow. 


planes of the mixing zone respectively, as shown in 
Fig. 2. If the pressure at the exit of the isentropic region, 
p c is specified then, given the known stagnation there and 

A 

the area ratio of the exit to the tube, — - the conditions at 

At 

the exit plane are completely known. Equations 8-10 may 
then be used to find p 0 , po, and uq. 

The conditions pi, p h and ui in the first computational 
cell of the tube are also known. The difference between p 0 
and pi will give rise to either a compression or expansion 
wave travelling to the right. The pressure ratio across the 
wave and the conditions in the first cell are sufficient 
information to analytically determine the velocity behind 
it, Uq* using Riemann invariants 8 . Since nothing has been 
said regarding the choice of p e , it is not expected that Uo* 
will be the same as Uo calculated using equations 8-10. A 
function may be defined however, as 

y(p e )=uo-uo* (11) 
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The proper choice of p e is that for which y is identically 
zero. This cannot be found analytically but may be 
obtained using a convenient numerical root finding 
technique (e.g. the false point method 9 ). The values of 
p 0 , Po, and uo obtained through the solution of equation 
1 1 are then assigned to the left image cell of the 
computing space. 

For outflow, consider Fig. 2(b). Here, the external region 
is modeled using a constant, known static pressure, p e . 
Adjacent to this is a fictitious zone in which the flow is 
isentropic and which, like the mixing zone above, adjusts 
instantaneously to flow conditions. Again, conditions in 
the first cell of the computing space are known. Under 
these modeling assumptions, for any value of po chosen a 
value of uo may be found to satisfy characteristic relations 
between the first cell and the image cell. Similarly, a 
value of uo* may be found to satisfy steady isentropic 
relations between the exit plane and the image cell. 
Through iteration, in a manner similar to the inflow 
scenario, a value of p 0 may ultimately be obtained for 
which uo and uq are identical. Detailed descriptions of 
similar boundary conditions may be found in Ref. 7. 

Choked Flow 

In this type of boundary condition, the quantity p 0 u 0 is 
specified along with the external stagnation temperature 
T 0 . The unknown quantities p 0 and uo required for the 
image cell are found through Riemann invariants and 
iteration. 

Wall 

For this boundary condition the quantities po, Po. and u 0 
are simply assigned the values p*, pi and -uj 
respectively. 


sector 2 




Figure 3. — Cross-sector compatibility schematic. 


Crossing Sectors/Area Change 

Compatibility between sectors is ensured by combining 
a fully-open type boundary condition for one sector and 
a partially-open boundary condition of another. For 
example, consider the geometry shown in Fig. 3 where 
the flow is from left to right. A guess is made for the 
exit pressure of sector 1. The state of the fully-open 
image cell for this sector is used to calculate stagnation 
properties. These, in turn, are used to determine the 
image cell state for the partially-open configuration of 
sector 2. The area of each sector is known. Thus, for the 
assumed exit pressure of sector 1, the mass flux (puA) 
in the image cell of each sector is also known. The 
proper choice of sector- 1 exit pressure then is that 
which yields the same mass flux in both sectors. This 
must be found through iteration. 

If the geometry is and flow conditions are such that flow 
is from a large to a small area sector, a similar iteration 
can be established using the fully-open inflow and 
partially-open outflow boundary conditions described 
above. 

At present, the area change aspect of the simulation is 
limited to unchoked flow; however, the method just 
described can easily be modified to accommodate 
choked flow with abrupt area change. 

Acoustic Validation 

As a check on the validity of the sectored-one- 
dimensional approach just described, a non-reactive, 
purely Euler calculation (i.e. the source vector 
identically zero) was performed on a tube with the 
geometry shown as the dashed line in Fig. 4. This 
geometry approximates the cross sectional area of a 



Figure 4. — Sectored and Q-l-D approximations to 
the area profile of a conical section. 
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conical section shown as a solid line. Also shown in the 
figure, as a dotted line, is the profile that was used for a 
quasi-one-dimensional (Q-l-D) calculation. Details of 
the Q-l-D formulation are not presented here. The 
calculation results are included for comparative 
purposes. The fundamental non-dimensional frequency 

for this shape ( f ' = ■—■ ) is known to be f=0.5 with 
a 

harmonics of 1.0, 1.5, etc 9 . Note that this is in contrast 
to a tube of uniform cross section and closed at one end, 
where the fundamental and harmonic frequencies are 
0.25, 0.75, 1.25, etc. If the conical tube simulations are 
excited by a pressure impulse at the open ends from a 
state of uniform pressure, temperature and zero 
velocity, the responses should contain the fundamental 
mode and several harmonics. The results of just such an 
excitation are shown in the Power Spectrum of Fig. 5. 
The computations were done with 100 numerical cells. 
The non-dimensional duration of both Q-l-D and 
Sectored simulations was 50 time units. The time 
interval for the data used in power spectrum analysis 
was 0. 1 time units. The excitation was a pressure spike 
at the open end 10% above ambient pressure, lasting 
0.5 Ume units. Pressure from the first interior cell, 
measured from the left end, was used to make the 
spectrum. The computed results are very similar for the 
two methods. Both Q-l-D and Sectored spectra are 
correctly dominated by the modes corresponding to the 
conical profile that was being simulated. While this not 
conclusive validation of the method just presented, it 
does provide confidence in the results that follow. 

Pulse Jet Validation 

As a second validation demonstration of the simulation 
method a Pulse Jet was considered. These propulsion 
devices, which capitalize on the coupling between 



0 0.5 1 1.5 2 

Frequency 


Figure 5. — Computed sectored and Q-l-D power 
spectral densities of pressure in the first cell of an 
approximated conical area profile using a single 
0.5 time unit exit pressure spike excitation. 


pressure fluctuation and heat-release, have been around 
for some time. They are generally very inefficient but 
also very simple in that they have no moving parts 
whatsoever. Commercial devices have been 

constructed. One such device, built by the SNECMA 
Company of France, has been documented in the 
literature 11 . The geometry is critical to performance; 
however, it should be possible to obtain large sustained 
oscillations with a geometry that is similar to a pulse 
combustor (these devices typically operate with 
pressure oscillations well above normal acoustic levels). 
Such a simplified geometrical layout, appropriate for 
the S-l-D code is shown in Fig. 6. Numbers along the 
top of this figure represent cross sectional areas scaled 
by the combustion section. The numbers along the 
bottom represent fractions of the total length. 

This simulation was run until limit cycle behavior was 
achieved. Limit cycle behavior is defined here as a 
repeating cycle with no change in peak-to-peak 
amplitude. The resulting wave pattern shown in Fig. 7. 
The figure shows contour shades of non-dimensional 
pressure temperature, Mach Number, and rate of heat 
release. Reference conditions for the non- 

dimensionalization correspond to Sea-Level- Standard. 
The horizontal direction corresponds to distance along 
the pulse jet. The vertical direction corresponds to non- 
dimensional time. Numbers to the left of each contour 
represent the maximum and minimum value found in the 
x-t space. It is noted that the simulation had run for at 
least one 200 time units before this figure was made. In 
this figure, the total time is 5.0 non-dimensional units. It 
is clear that strong, self-sustained oscillations are present 
as expected. Peak-to-peak pressure oscillations of 52% of 
the mean value were computed in the combustion 
chamber. The frequency of oscillations corresponds to 
that reported in the literature. The simulation was made 
with 200 cells. The mean fuel/air ratio for these devices is 
unknown; however, they are designed to run under 
standard atmospheric inlet conditions. It was assumed 
that the peak temperature allowed is 2300 R. The reactant 
source term, m r in the source vector (Eqn. 5) was 
applied to a single numerical cell at the location shown in 
Fig. 6 and was adjusted to achieve this value. 
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Figure 6. — Sectored approximation of a 
SNECMA pulse jet. 
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Figure 7. — Contours of non-dimensional pressure, temperature, 
Mach number, and heat release rate during several cycles 
of a SNECMA pulse jet simulation. 


As expected, this and subsequent simulation results to 
be presented satisfied Rayleigh's Criteria 12 for 
instability. Figure 8 shows the integral 

i 

Ra = q 0 jRp'dx (12) 

o 

for the Pulse Jet simulation over the approximately 
three cycle time period illustrated in Fig. 7. The primes 
in Eqn. 12 refer to fluctuations from the time-mean 
value of reaction rate and pressure. Rayleigh’s Criteria, 
which is essentially a measure of the phase relationship 
between heat-release and pressure fluctuation, requires 
that the integral of Ra over one oscillation be positive, 
which it clearly is. 

It is interesting to note in passing that Pulse Jets are 
decidedly non-premix devices, yet were successfully 
simulated here. This may suggest utility of the method 
just presented beyond the original premix restriction. 



Figure 8. — Computed Rayleigh Index for the 
SNECMA pulse jet simulation. 


Combustor Simulation Results 
In a further effort to demonstrate the capabilities 
afforded by the method described above, a simulation of 
an experimental combustion rig exhibiting known 
instabilities was implemented. The instabilities 
exhibited by the rig were associated with sections 
beyond the expected injector and combustion zone. As 
such it was necessary to include the downstream plenum 
and cooling ring in the computing domain. The 
geometry simulated is shown in Fig. 9. For the results to 
be shown, 348 numerical cells were used. 



x (inches) 

Figure 9. — Combustion rig geometry. 
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The left side boundary conditions (x=0.0 in.) were 
modeled as fully-open, constant pressure. The right side 
boundary conditions (x=99.6 in.) were those for a solid 
wall. Flow exited the rig through a tee-section pipe 17 
in. from the right end. This was modeled in the 
simulation using an orifice type equation representing 
m s in the source vector, distributed over a number of 

computational cells. Flow through each orifice is driven 
by the ratio of pressure in the numerical cell to that of 
an imposed back pressure. 


Table 2. — Combustor simulation parameters. 


£ t /Re* 

0.00056 

Kq 

26.0 

flame 

27.6 in. 

Sc t 

1.0 

Pr, 

1.0 

Ti 

1.25 

qo 

7.0 

Pback 

.54932 

At 

.00133 


The cooling spray was modeled using the simple 
conductive Eqn. 7 applied to all computational cells in 
the domain beyond (to the right of) x=59.62 in. The value 
used for T^ in this equation was the measured exhaust 
temperature of approximately 760 R. The value of a was 
adjusted until the gas temperature at the beginning of the 
tee-section was equal to the measured exhaust 
temperature. This type of simplistic model implies that 
the dominant effect of the cooling spray is the changing 
of the gas temperature (e.g. speed of sound). This may 
not always be the case; however, the model seems to be 
effective for the rig results presented here. 

Like the Pulse Jet simulation, fuel injection was 
modeled with the m r term of the source vector applied 

to a single numerical cell at the location shown in 
Fig. 9. The source term was held at a constant value for 
each operating point to be described. This implies a 
constant fuel flow rate, which is appropriate due to the 
large fuel supply pressure of the rig. 

Two operating points were simulated that exhibited two 
distinct frequency instabilities. The conditions for the 
operating points are listed in Table 1 . The subscript 3 in 
this table represents an upstream location near the inlet. 
The subscript 4 represents a downstream location near 
x=40.0 in. It is noted that conditions at Operating Point 
2 also represent the reference conditions used for non- 
dimensionalization. Thus, any non-dimensional results 
presented may be dimensionalized using these values. 

Other relevant parameters of the simulations are listed in 
Table 2. The first three values in this table were chosen, ad- 
hoc, because they produced instability in the simulation. 
The turbulent Schmidt and Prandtl numbers, and the 


Table 1. — Combustor operating points. 


Pnt. 

P 3 

Psia 

T 3 

R 

m 

lbm/s 

Y 

A/F 

t 4 

R 

1 

155 

1101 

0.665 

1.35 

0.038 

3307 

2 

244 

1361 

L34 

1.35 

0.024 

2858 


ignition temperature were chosen based on engineering 
estimates. The value of q 0 was somewhat arbitrary. The 
value used, in conjunction with m r , determines the 

available chemical energy. The time step, At, was chosen 
based on numerical stability considerations. 

The simulated position of the flame (reaction front) is 
determined by a prescribed distribution of the various 
diffusion coefficients (species, temperature, and 
momentum) as described earlier. The existence of an 
unstable operating point is often quite sensitive to the 
flame position both in terms of its proximity to the 
injection point, and in terms of its overall location in the 
combustor. As such, if an instability is observed, 
moving the flame position slightly restores steady-state 
operation. This technique was used in the simulation to 
achieve initially stable operation at the conditions listed 
in Table 1. 


For reference, stable distributions of temperature and 
Mach Number for Operating Point 1 are shown in 
Fig. 10. Pressure and non-dimensional reaction rate 
(R in Eqn. 6) are shown in Fig. 11. The slightly wavey 
appearance of the Mach Number and pressure at the left 
side of the figures is due to the deliberate introduction 
of small random pressure oscillations at the inlet. These 
had a maximum value of ±0.5% of the mean. They 
were intended to simulate noise from unknown sources 
in the system. 



Figure 10. — Steady-state distributions of 
temperature, and Mach number at 
operating point 1. 
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x (inches) 


Figure 11. — Steady-state distributions of pressure, 
and non-dimensional reaction rate at 
operating point 1. 


It should be noted that the rig geometry of the two 
operating points was slightly different. For Operating 
Point 1, the narrow region where the fuel is injected was 
1.0 in. shorter than is shown in Fig. 9. The large cross 
section sector preceding this was 1.0 in. longer than is 
shown. Thus the overall length was the same. For 
Operating Point 2, the geometry was as shown in Fig. 9. 
These slight differences were included in each 
simulation. 

Operating Point 1 

With the flame position properly placed, self-excited, 
unstable operation commenced, and eventually reached 
limit cycle behavior. Limit cycle behavior is defined 
here as a repeating cycle that exists after at least 4 
seconds of simulated time with no significant change in 
peak-to-peak amplitude. The computed peak-to-peak 
amplitude of the pressure fluctuations at x=42.1 in. (the 
location of a pressure transducer) was approximately 
7.7% of the mean. This is close to the 8% value 
measured experimentally. Contour plots of non- 
dimensional pressure and velocity fluctuations 
(deviations from the time-mean values) are presented in 
Fig. 12 for approximately 2 oscillatory cycles. These 
are presented in order to show the mode shape of the 
instability. It would appear that this limit cycle is 
primarily the quarter wave mode of the combustor 
section proper; however, other modes are present. 

Measured and computed Power Spectral Density 
distributions (PSD) of pressure at the experimental 
transducer location of x=42.1 in. are presented in 
Fig. 13. The signal used for the computed results was of 
0.116 seconds duration at a 10 khz. sampling rate. For 
the measured results the signal duration was at least 
10 seconds at the same 10 khz. sampling rate. The 
match between measured and computed power spectra 
is, in many ways, remarkable. The critical similarities, 
of course, are the frequency and power density 
associated with the instability. For reference, the 
measured dominant frequency was 275 hz. The 


Pressure 



o.o 0.5 1.0 

x/L 


Velocity 



0.0 0.5 1.0 

x/L 


Figure 12. — Contours off non-dimensional pressure 
and velocity fluctuations for 2 cycles of the 
operating point 1 limit cycle. 



Figure 13. — Power spectral density distribution of 
pressure for x=42.1 inches during the operating 
point 1 limit cycle. 

computed value was 291 hz. Although differing 
significantly in magnitude, the presence of higher 
frequency modes can be seen at approximately 570 and 
850 hz. in both the computed and measured spectra. 
This suggests that the subtle acoustic properties of the 
rig are reasonably reproduced in the simulation. It is 
interesting to note that the 850 hz. signal corresponds to 
the 3/4 wave mode of the combustor section proper. 
The 550 hz. peak however seems to be unrelated. 
Another striking similarity between measured and 
computed power spectra is the visible low frequency 
peak of approximately 188 hz. What makes this striking 
is that first, it is the dominant mode found at Operating 
Point 2 and second, that it is not a subharmonic of the 
275 hz. mode. 

Beyond 1000 hz., there are large disparities between 
measured and computed power spectra. The reasons for 
this are unclear. One possible explanation is that there 
are other acoustic sources in the vicinity of the 
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experimental pressure transducer such as a cooling 
water spray bar, and vorticies shed from recirculation 
zones that are not modeled in the simulation. 

Operating Point 2 

As with Operating Point 1, the simulation successfully 
reproduced the observed instability at approximately the 
correct frequency and amplitude. The measured peak- 
to-peak pressure oscillations were 3.6 % of the mean. 
The computed oscillations were 3.1 % of the mean. The 
measured dominant frequency was 175 hz. while the 
computed value was 188 hz. Figure 14 shows the same 
type contour plot as Fig. 12 in order to illustrate the 
dominant oscillatory mode. It is clear that this mode is 
completely different in form from that of Fig. 12 and 
that the downstream plenum of the rig is acoustically 
active. This mode is the main reason that it was 
necessary to simulate such a large portion of the rig. 

The measured and computed power spectra for this 
operating point are shown in Fig. 15. As before, the 
computed and measured results agree well with respect 
to amplitude and frequency of the dominant instability. 
The qualitative agreement is also quite good in some 
respects. In particular, both show the dominant low 
frequency mode, and the virtual absence of the 275 hz. 
mode found at Operating Point 1. Both also show a 
second peak near 380 hz. The source of this is 
unknown. The measured results show a peak at 450 hz. 
which is not seen in the computed results. The 
explanation for this is unclear at the present time. 

Discussion 

The results of the simulation just presented are 
encouraging in as much as they successfully duplicated 
self-excited instabilities observed in a rig with relatively 
complex geometry. It is also encouraging that a 
simulation of this combustor rig lasting 0.116 seconds 
of computed time took 188 seconds of CPU time 
running on a Sun Ultra 2 Workstation. 


Pressure 



0.0 0.5 1.0 

x/L 
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Figure 14. — Contours off non-dimensional pressure 
and velocity fluctuations for 2 cycles of the 
operating point 2 limit cycle. 


10 ° 



Frequency (hz.) 


Figure 15. — Power spectral density distribution of 
pressure for x=42.1 inches during the operating 
point 2 limit cycle. 

delays rather than to pressure induced variations in the 
reaction rate itself. The latter mechanism has been 
presented before in purely one-dimensional combustor 
simulations 3 . The former mechanism is a more common 
explanation for combustors of the lean-premix type 
described here. 


It is interesting to note, with regard to the objective of 
successfully capturing the physical phenomena of the 
instability, that the use of a fuel injection model was 
necessary to obtain the rig results presented. It is 
possible to inject the fuel (reactant) at the boundary 
(inlet) of the simulation rather than at an interior cell. 
This is done by simply specifying the reactant fraction 
at the boundary. When this was done for the rig 
geometry just described, maintaining exactly the same 
flow rates and temperatures, the instability disappeared 
completely (this was not the case for the Pulse Jet 
simulation). This suggests that the coupling mechanism 
between heat release and pressure is predominantly 
related to variations of mixture ratio and convection 


Limitations 

Restrictions regarding the applicability of the simulation 
method just described were mentioned in the 
introduction (premix combustors with longitudinal 
modes only). However, it is important to note several 
other important limitations. 

Diffusivitv Distribution 

The particular rig simulation was, as mentioned, very 
sensitive to the prescribed distribution and magnitude of 
diffusivity. Recall that this distribution essentially 
determines the location of the reaction zone. Given the 
apparent dependence of the instabilities on convection 
time between the fuel injection point and the reaction 
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zone, the sensitivity is not surprising. The distribution 
however, is completely heuristic at the time of this 
paper and as such, it must be concluded that it was 
'tuned' in order to obtain results that were observed. In 
other words, had a slightly different distribution been 
chosen, no self-excited instability would have been 
observed. For example, if the flame position for 
Operating Point 1 is moved downstream just 0.3 in. 
(one grid point), the 292 hz. mode disappears and the 
188 hz. mode dominates. If it is moved down 0.6 in. 
instabilities of any kind vanish. Furthermore, when the 
simulation was run at a third operating point, nearly 
identical to Operating Point 1 but with a higher inlet 
temperature, it produced essentially the same 292 hz. 
mode while experimentally no instability was detected. 
To some degree, this sensitivity may be a result of the 
particular rig simulated. Until a more physically based 
model can be developed for diffusivity however, it 
would seem that the use of this simulation as a reliable 
predictor of instability is premature. This is particularly 
true since, as stated, the use of diffusivity is itself a 
model of the effects of recirculating flow. 

Another difficulty with the diffusivity distribution is that 
it cannot cross sectors. Diffusion terms in the governing 
equations can only be applied to interior numerical 
cells. Thus, since each sector is separated by a 
boundary, diffusion processes must begin and end 
within a sector. 

Flame Structure 

Related to the flame position or reaction zone, it should 
be kept in mind that, despite the geometric 
simplifications described herein, the actual flame 
structure in most lean, premix combustors is highly 
three dimensional. It does not closely approach the 
planar form assumed in this simulation. It is possible 
however, that an integrated average of the muti- 
dimensional combustion/recirculation process would 
closely match that of the present simulation. This should 
be assessed using high fidelity numerical codes. 

Multiple area Changes 

Because the boundary condition routines associated 
with the sectored approach described here are highly 
iterative, their computational overhead is large. As 
such, the method becomes prohibitively time consuming 
if the number of area changes in a simulation is large. 

Conclusions 

The sectored-one-dimensional, numerical, reacting flow 
solver technique described in this paper reliably 
simulates the behavior of thermo-acoustic instabilities 
observed in premix combustors with complex cross- 
sectional geometries. The fact that it is relatively simple 


(e.g. one-dimensional), yet reasonably accurate make it 
useful for characterization and testing of control 
strategies. Furthermore, modeling and implementation 
of control actuation, such as fuel flow perturbations, are 
straightforward in the framework of the simulation. 
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Appendix: Difficulties of 0-1 -D Formulations 


It would appear at first that a Quasi-One-Dimensional 
numerical approach would be practical for the 
combustor geometries described in this paper. Indeed, a 
reactive Q-l-D code was developed by the author as an 
initial approach to modeling the experimental rig 
described herein and work continues on this approach. 
However, there are difficulties presented by this 
approach which, to date, the author has not been able to 
successfully resolve. All of the difficulties revolve 
around the implementation of the momentum source 

dA i « . 

term p — . The appropnate manner in which to 
dx 

evaluate this term numerically is somewhat ambiguous 
when the area gradients are large. Among other things, 
it seems to depend on the numerical integration scheme 
used. Furthermore, even if it is evaluated in a manner 
that produces a stable scheme, other difficulties may 
arise. These include, among others: 

1. The need for many grid points in vicinities of area 
change. 

2. Non-physical behavior such as local entropy loss 
(gains in total pressure) and/or unpredictable total 
pressure losses related to the area change. 

3. Changing the acoustic properties of the system due 
to approximating abrupt area changes with something 
more gradual. 


Additionally, Q-l-D formulations often result in flame 
locations precisely in vicinities of large area gradients. 
Since thermo-acoustic instabilities result from coupling 
between heat release and pressure perturbations, it may 
become unclear whether simulated instabilities are 
physically real or the result of numerics. 

In any case, as mentioned in an earlier footnote, a Q-l- 
D simulation has not yet been implemented that 
successfully simulates the rig instabilities just described. 
They have however, been used to successfully simulate 
some Pulse Jet engines. 
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